#include "package.h"

double Package::cal(State_TP &TP_f, State_TP &TP_i) {
  double result = 0.0;
  double part1, part2, part3, part4;
  if (TP_f.J_2 != TP_i.J_2 || TP_f.T_2 != TP_i.T_2)
    return result;
  int ch_TP = TP.FindC(TP_f.J_2, TP_f.T_2);
  int ch_NL = jacNL.FindC(TP_f.J_2, TP_f.T_2);
  int Ep = TP_f.E;
  int TPp_index = TP_f.index - TP.channel[ch_TP].indexE_beg[Ep];
  int E = TP_i.E;
  int TP_index = TP_i.index - TP.channel[ch_TP].indexE_beg[E];

  int i_beg = jacNL.channel[ch_NL].indexE_beg[Ep];
  int i_end = jacNL.channel[ch_NL].indexE[Ep];
  int j_beg = jacNL.channel[ch_NL].indexE_beg[E];
  int j_end = jacNL.channel[ch_NL].indexE[E];
  // cout<<endl;
  // cout<<"channel_NL : "<<ch_NL<<" Ep : "<<Ep<<" E : "<<E<<endl;
  // cout<<"i_beg : "<<i_beg<<" i_end : "<<i_end<<" j_beg : "<<j_beg<<" j_end :
  // "<<j_end<<endl;

  part1 = 0.0;
  for (int i = i_beg; i <= i_end; i++) {
    part2 = 0.0;
    for (int j = j_beg; j <= j_end; j++) {
      if (jacNL.state_NL[ch_NL][i].N_0 != jacNL.state_NL[ch_NL][j].N_0 ||
          jacNL.state_NL[ch_NL][i].L_0 != jacNL.state_NL[ch_NL][j].L_0)
        continue;
      int i_index = i - i_beg;
      int j_index = j - j_beg;
      // cout<<i_index<<"\t"<<j_index<<endl;
      if (trans.Tval_vec_flag[ch_TP][Ep][i_index][TPp_index] != 1 ||
          trans.Tval_vec_flag[ch_TP][E][j_index][TP_index] != 1) {
        cout << "Wrong @ Package::cal2" << endl;
        cout << "trans.Tval_vec_flag[ch_TP][Ep][i_index][TPp_index] != 1"
             << endl;
        exit(0);
      }
      double Tp = trans.Tval_vec[ch_TP][Ep][i_index][TPp_index];
      double T = trans.Tval_vec[ch_TP][E][j_index][TP_index];
      double TpT = Tp * T;
      if (TpT == 0.0)
        continue;

      int ii = jacNL.state_NL[ch_NL][i].index;
      int jj = jacNL.state_NL[ch_NL][j].index;
      int E_ii = jacNL.state_NL[ch_NL][i].E;
      int E_jj = jacNL.state_NL[ch_NL][j].E;
      int ch_jac_ii = jacNL.state_NL[ch_NL][i].JT_channel;
      if (ch_jac_ii != jacNL.state_NL[ch_NL][j].JT_channel)
        continue;
      int ch_jac = ch_jac_ii;
      int x_beg = jac.channel[ch_jac].indexE_beg[E_ii];
      int x_end = jac.channel[ch_jac].indexE[E_ii];
      int y_beg = jac.channel[ch_jac].indexE_beg[E_jj];
      int y_end = jac.channel[ch_jac].indexE[E_jj];
      int ii_index = ii - x_beg;
      int jj_index = jj - y_beg;
      part3 = 0.0;
      for (int x = x_beg; x <= x_end; x++) {
        part4 = 0.0;
        for (int y = y_beg; y <= y_end; y++) {
          int x_index = x - x_beg;
          int y_index = y - y_beg;
          double Dpx = anti.vec_D[ch_jac][E_ii][ii_index][x_index];
          double Dy = anti.vec_D[ch_jac][E_jj][jj_index][y_index];
          double cont_term = cont.vec[ch_jac][x][y];
          double ope_term = ope.vec[ch_jac][x][y];
          double part4_t = Dpx * Dy * (cont_term + ope_term);
          part4 += part4_t;
        }
        part3 += part4;
      }
      part2 += part3 * TpT;
    }
    part1 += part2;
  }
  result = 3 * part1 * 6;
  // cout<<result<<endl;
  return result;
}

void Package::build(int TP_ch) {
  int size = TP.state[TP_ch].size();
  vec.resize(1);
  vec[0].resize(size);

  int i = 51;
  int j = 32;
  double result = cal2(TP.state[TP_ch][i], TP.state[TP_ch][j]);
  // cout << i << "\t" << j << "\t" << result << endl;
  /*
  for (int i = 0; i < size; i++) {
    // cout<<"i : "<<i<<endl;
    vec[0][i].resize(size);
    for (int j = 0; j < size; j++) {
      double result = cal2(TP.state[TP_ch][i], TP.state[TP_ch][j]);
      vec[0][i][j] = result;
      // if(result != 0.0)
      // cout<<i<<"\t"<<j<<"\t"<<result<<endl;
    }
  }*/
}

void Package::build() {
  LEC_CONST LECs;
  this->LEC_CD_P = LECs.LEC_CD;
  this->LEC_CE_P = LECs.LEC_CE;
  this->LEC_C1_P = LECs.LEC_C1;
  this->LEC_C3_P = LECs.LEC_C3;
  this->LEC_C4_P = LECs.LEC_C4;
  this->LEC_gA_P = LECs.LEC_gA;

  int ch_size = TP.channel.size();
  vec.resize(ch_size);
  vec_flag.resize(ch_size);
  Pgs_Bar pgs;
  pgs.init(ch_size, 30);
  this->vec_S.init_ch(ch_size, 1);
  this->vec_flag_S.init_ch(ch_size, 0);
  cout << "Begin to build Package : " << endl;
  for (int TP_ch = 0; TP_ch < ch_size; TP_ch++) {
    pgs.PrintBar(TP_ch);

    int size = TP.state[TP_ch].size();
    vec[TP_ch].resize(size);
    vec_flag[TP_ch].resize(size);
    // cout<<"channel : "<<TP_ch<<endl;
    //#pragma omp parallel for schedule(dynamic, 1)
    for (int i = 0; i < size; i++) {
      vec[TP_ch][i].resize(size);
      vec_flag[TP_ch][i].resize(size);

      for (int j = 0; j < size; j++) {
        vec[TP_ch][i][j] = 0;
        vec_flag[TP_ch][i][j] = 0;
      }
    }
    // int flag_i=-1;
    // if(i>0)flag_i =
    // TP.TzDiff(TP.state[TP_ch][i],TP.state[TP_ch][i-1]);//flag =1;
    // degenerate of Tz
    // if(flag_i == 1 ){
#pragma omp parallel for schedule(dynamic, 1)
    for (int k = 0; k < TP.TPconf[TP_ch].size(); k++) {
      int i = TP.TPconf[TP_ch][k].F;
      int j = TP.TPconf[TP_ch][k].I;
      /*for (int j = 0; j < size; j++) {
        int flag = -1;

        if (TP.state[TP_ch][i].F > 0 && TP.state[TP_ch][j].I > 0)
          flag = 1;
        // if (TP.state[TP_ch][y].F > 0 && TP.state[TP_ch][x].I > 0)
        //  flag = 1;
        if (flag < 0)
          continue; */
      // vec[TP_ch][i][j]=vec[TP_ch][i-1][j];
      //}
      // }else{
      // for(int j=0; j<size; j++){
      // int flag_j=-1;
      // if(j>0)flag_j =
      // TP.TzDiff(TP.state[TP_ch][j],TP.state[TP_ch][j-1]);//flag =1;
      // degenerate of Tz
      // if(flag_j == 1){
      // vec[TP_ch][i][j]=vec[TP_ch][i][j-1];
      //}else{
      // cout << "TP_ch : " << TP_ch << "\t i = " << i << "\t j = " << j <<
      // endl;
      double result = this->cal2(TP.state[TP_ch][i], TP.state[TP_ch][j]);
      vec[TP_ch][i][j] = result;
      vec_flag[TP_ch][i][j] = 1;
      // if(result != 0.0)
      // cout<<i<<"\t"<<j<<"\t"<<result<<endl;
      //}
    }
    this->vec_S.build_ch(this->vec[TP_ch], TP_ch);
    this->vec[TP_ch].clear();
    this->vec_flag_S.build_ch(this->vec_flag[TP_ch], TP_ch);
    this->vec_flag[TP_ch].clear();
  }
  pgs.EndBar();
  cout << "End to build Package : " << endl;
  this->vec.clear();
  this->vec_flag.clear();
  // cout << " Start to build vec_S.build_ch " << endl;
  // this->vec_S.build_ch(this->vec);
  // this->vec_flag_S.build_ch(this->vec_flag);
  // cout << "End to build Package vec_S " << endl;
}

void Package::build_JacAvec() {
  cout << " Begin build_JacAvec " << endl;
  int size_ch = jac.state.size();
  // size_ch=1;
  JacA_vec.resize(size_ch);
  JacA_vec_flag.resize(size_ch);

  for (int ch = 0; ch < size_ch; ch++) {
    int size = jac.state[ch].size();
    // size=2;
    JacA_vec[ch].resize(size);
    JacA_vec_flag[ch].resize(size);
    //#pragma omp parallel for schedule(dynamic, 1)
    for (int i = 0; i < size; i++) {
      JacA_vec[ch][i].resize(size);
      JacA_vec_flag[ch][i].resize(size);
      for (int j = 0; j < size; j++) {
        JacA_vec[ch][i][j] = 0;
        JacA_vec_flag[ch][i][j] = -1;
      }
    }
#pragma omp parallel for schedule(dynamic, 1)
    for (int k = 0; k < jac.Jacconf[ch].size(); k++) {
      int i = jac.Jacconf[ch][k].F;
      int j = jac.Jacconf[ch][k].I;
      // for (int j = 0; j < size; j++) {
      //   if (jac.state[ch][i].F < 0 || jac.state[ch][j].I < 0)
      //     continue;
      int E_i = jac.state[ch][i].E;
      int E_j = jac.state[ch][j].E;
      int x_beg = jac.channel[ch].indexE_beg[E_i];
      int x_end = jac.channel[ch].indexE[E_i];
      int y_beg = jac.channel[ch].indexE_beg[E_j];
      int y_end = jac.channel[ch].indexE[E_j];
      double result = 0.0;
      for (int x = x_beg; x <= x_end; x++) {
        double D_ix = anti.vec_D[ch][E_i][i - x_beg][x - x_beg];
        if (D_ix == 0.0)
          continue;
        for (int y = y_beg; y <= y_end; y++) {
          double D_jy = anti.vec_D[ch][E_j][j - y_beg][y - y_beg];
          if (D_jy == 0.0)
            continue;
          // cout<<x<<"\t"<<y<<endl;
          // cout<<D_ix<<"\t"<<D_jy<<"\t"<<cont.vec[ch][x][y]<<endl;
          // !!!! Be careful here !!!
          /*double cont_val = cont.quick(ch, x, y);
          double ope_val = ope.quick(ch, x, y);
          double tpe_val = tpe.quick(ch, x, y);
          double mat = cont_val + ope_val + tpe_val; */
          if (tpe.vec_flag_tot[ch][x][y] != 1 &&
              (jac.state[ch][x].Par == jac.state[ch][y].Par)) {
            cout << "Wrong happened at Package::build_JacAvec 1" << endl;
            cout << "tpe.vec_flag_tot[ch][x][y]!=1" << endl;
            cout << "channel : " << ch << "\t x : " << x << "\t y : " << y
                 << "\t tpe.vec = " << tpe.vec[ch][x][y] << endl;
            exit(0);
          }
          double mat =
              cont.vec[ch][x][y] + ope.vec[ch][x][y] + tpe.vec[ch][x][y];
          // double mat = tpe.vec[ch][x][y] + ope.vec[ch][x][y];
          // double mat = tpe.vec[ch][x][y];
          result += D_ix * D_jy * mat;
        }
      }
      JacA_vec[ch][i][j] = 6 * result;
      JacA_vec_flag[ch][i][j] = 1;
      // cout<<JacA_vec[ch][i][j]<<endl;
      //}
    }
  }
  cout << " End build_JacAvec " << endl;
}

void Package::print_JacAvec(std::string filename) {
  std::ofstream file;
  file.open(filename);
  file << "!!! Cation the print is 3 times of vec result" << endl;
  int size_ch = JacA_vec.size();
  for (int ch = 0; ch < size_ch; ch++) {
    int size = JacA_vec[ch].size();
    file << " Channel : " << ch << endl;
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
        if (JacA_vec[ch][i][j] != 0)
          file << i << "\t" << j << "\t" << JacA_vec[ch][i][j] * 3 << endl;
      }
    }
  }
  file.close();
}

double Package::cal2(State_TP &TP_f, State_TP &TP_i) {
  // use JacA_vec to calculate
  double result = 0.0;
  if (TP_f.J_2 != TP_i.J_2 || TP_f.T_2 != TP_i.T_2)
    return result;
  int Par_f = phase_cal(TP_f.E);
  int Par_i = phase_cal(TP_i.E);
  if (Par_f != Par_i)
    return result;

  // if (TP_f.T_2z != TP_i.T_2z)
  //  return result;
  // if(TP_f.Tab_2 != TP_i.Tab_2)return result;
  // cout << "good ? 1" << endl;
  int ch_TP = TP.FindC(TP_f.J_2, TP_f.T_2);
  int ch_NL = jacNL.FindC(TP_f.J_2, TP_f.T_2);
  int Ep = TP_f.E;
  int TPp_index = TP_f.index - TP.channel[ch_TP].indexE_beg[Ep];

  int E = TP_i.E;
  int TP_index = TP_i.index - TP.channel[ch_TP].indexE_beg[E];

  // cout<<"good?"<<endl;

  int i_beg = jacNL.channel[ch_NL].indexE_beg[Ep];
  int i_end = jacNL.channel[ch_NL].indexE[Ep];
  int j_beg = jacNL.channel[ch_NL].indexE_beg[E];
  int j_end = jacNL.channel[ch_NL].indexE[E];
  if (std::abs(i_end - i_beg) < 1 || std::abs(j_end - j_beg) < 1) {
    return result;
  }
  /*cout << endl;
  cout << " ====  E_p = " << Ep << "\t E = " << E << " ==== " << endl;
  cout << "\t i_beg = " << i_beg << "\t i_end = " << i_end
       << "\t\t j_beg = " << j_beg << "\t j_end = " << j_end << endl;
 */
  // cout << "good ? 2" << endl;
  // double part1,part2,part3,part4;
  double part1 = 0.0;
  for (int i = i_beg; i <= i_end; i++) {
    int i_index = i - i_beg;
    if (trans.Tval_vec_flag[ch_TP][Ep][i_index][TPp_index] != 1) {
      cout << "Wrong @ Package::cal2" << endl;
      cout << "trans.Tval_vec_flag[ch_TP][Ep][i_index][TPp_index] != 1" << endl;
      exit(0);
    }
    double T_p = trans.Tval_vec[ch_TP][Ep][i_index][TPp_index];
    // cout << "\t good ? 2 - 1 - 1" << endl;
    int ii = jacNL.state_NL[ch_NL][i].index;
    int ch_jac_ii = jacNL.state_NL[ch_NL][i].JT_channel;
    // if(T_p ==0.0)continue;
    // cout << "\t good ? 2 - 2" << endl;
    double part2 = 0.0;
    for (int j = j_beg; j <= j_end; j++) {
      if (jacNL.state_NL[ch_NL][i].N_0 != jacNL.state_NL[ch_NL][j].N_0 ||
          jacNL.state_NL[ch_NL][i].L_0 != jacNL.state_NL[ch_NL][j].L_0)
        continue;

      int j_index = j - j_beg;
      // cout<<i_index<<"\t"<<j_index<<endl;
      if (trans.Tval_vec_flag[ch_TP][E][j_index][TP_index] != 1) {
        cout << "Wrong @ Package::cal2" << endl;
        cout << "trans.Tval_vec_flag[ch_TP][Ep][i_index][TPp_index] != 1"
             << endl;
        exit(0);
      }
      double T = trans.Tval_vec[ch_TP][E][j_index][TP_index];
      // if(T == 0.0)continue;
      int jj = jacNL.state_NL[ch_NL][j].index;
      if (ch_jac_ii != jacNL.state_NL[ch_NL][j].JT_channel)
        continue;
      /*
            if (std::abs(T * T_p) > -1.0) {
              cout << endl;
              cout << "-----  ch_jac_p = " << ch_jac_ii
                   << "\t ch_jac = " << jacNL.state_NL[ch_NL][j].JT_channel <<
         endl;
              cout << "\t i = " << i << "\t i_index = " << i_index
                   << "\t T_p = " << T_p << endl;
              cout << "\t j = " << j << "\t j_index = " << j_index << "\t T = "
         << T
                   << endl;
              cout << "\t *** ch_jac = " << ch_jac_ii << "\t ii = " << ii
                   << "\t jj = " << jj << "\t JacA = " <<
         JacA_vec[ch_jac_ii][ii][jj]
                   << endl;
            }                */
      // int ch_jac = ch_jac_ii;
      if (JacA_vec_flag[ch_jac_ii][ii][jj] != 1) {
        cout << "Wrong @ Package::cal2 JacA_vec_flag[ch_jac_ii][ii][jj] != 1 "
             << endl;
        cout << "ch_jac_ii = " << ch_jac_ii << "\t ii = " << ii
             << "\t jj = " << jj << endl;
        cout << "ch_NL = " << ch_NL << "\t i = " << i << "\t j = " << j << endl;
        exit(0);
      }
      part2 += T_p * T * JacA_vec[ch_jac_ii][ii][jj];
      // cout<<"part2 = "<<part2<<endl;
    }
    part1 += part2;
    // cout << "\t good ? 2 - 3" << endl;
  }
  result = 3 * part1;
  // cout << "good ? 3" << endl;
  // if(result != 0.0)cout<<result<<endl;
  return result;
}

void Package::print(std::string filename) {
  std::ofstream file;
  file.open(filename);
  int size_ch = vec_S.size_ch;
  for (int ch = 0; ch < size_ch; ch++) {
    int size = vec_S.row_num_ch[ch];
    file << "Channel : " << ch << endl;
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
        double val = vec_S.val(ch, i, j);
        // if (std::abs(vec[ch][i][j]) > 0.0001)
        if (std::abs(val) > 0.0001)
          file << i << "\t" << j << "\t" << val << endl;
      }
    }
  }
  file.close();
}

void Package::print_Tz(std::string filename) {
  std::ofstream file;
  file.open(filename);
  int size_ch = vec_Tz_S.size_ch;
  for (int ch = 0; ch < size_ch; ch++) {
    int size = vec_Tz_S.row_num_ch[ch];
    file << "Channel : " << ch << endl;
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
        double val = vec_Tz_S.val(ch, i, j);
        // if (std::abs(vec_Tz[ch][i][j]) > 0.000000001)
        if (std::abs(val) > 0.000001)
          file << i << "\t" << j << "\t" << val << endl;
      }
    }
  }
  file.close();
}

void Package::print_TzM(std::string filename) {
  std::ofstream file;
  file.open(filename);

  int size = vec_TzM.size();

  for (int i = 0; i < size; i++) {
    for (int j = 0; j < size; j++) {
      if (std::abs(vec_TzM[i][j]) > 0.0001)
        file << i << "\t" << j << "\t" << vec_TzM[i][j] << endl;
    }
  }
  file.close();
}

double Package::Check_Anti_T(SixJ_coeff &sixJ, SingleP &SP, int J2, int T2,
                             int flag, int index_f, int index_i) {
  double res = 0.0;
  int ch = TP.FindC(J2, T2);
  int a = TP.state[ch][index_i].a;
  int b = TP.state[ch][index_i].b;
  int c = TP.state[ch][index_i].c;
  int E = TP.state[ch][index_i].E;
  int T2z = TP.state[ch][index_i].T_2z;
  int J2ab = TP.state[ch][index_i].Jab_2;
  int T2ab = TP.state[ch][index_i].Tab_2;

  if (flag == 1) {
    // a <-> b
    int phase_t = SP.state[a].j_2 + SP.state[b].j_2 + J2ab + 2 + T2ab;
    int phase = 1 - 2 * ((phase_t / 2) % 2);
    res = vec[ch][index_f][index_i] * phase * -1;
    return res;
  } else {
    // b <-> c
    int size = TP.state[ch].size();
    double part1 = 0.0;
    int phase_t = SP.state[b].j_2 + SP.state[c].j_2 + J2ab + 2 + T2ab;
    int phase = 1 - 2 * ((phase_t / 2) % 2);
    double part1_t = (J2ab + 1) * (T2ab + 1);
    part1 = phase * std::sqrt(part1_t);
    double part2 = 0.0;
    for (int i = 0; i < size; i++) {
      double part2_t = 0.0;
      if (TP.state[ch][i].E != E || TP.state[ch][i].T_2z != T2z)
        continue;
      // if(TP.state[ch][i].E != E )continue;
      if (TP.state[ch][i].a != a)
        continue;
      if (TP.state[ch][i].b != c || TP.state[ch][i].c != b)
        continue;
      int J2ac = TP.state[ch][i].Jab_2;
      int T2ac = TP.state[ch][i].Tab_2;
      double part2_tt = (J2ac + 1) * (T2ac + 1);
      phase_t = J2ac + T2ac;
      phase = 1 - 2 * ((phase_t / 2) % 2);
      part2_t = phase * std::sqrt(part2_tt);
      part2_t = part2_t * sixJ.quick(SP.state[a].j_2, SP.state[b].j_2, J2ab, J2,
                                     SP.state[c].j_2, J2ac);
      part2_t = part2_t * sixJ.quick(1, 1, T2ab, T2, 1, T2ac);
      part2 += part2_t * vec[ch][index_f][i];
      // cout<<"\t"<<index_f<<"\t"<<i<<"\t"<<vec[ch][index_f][i]<<endl;
    }
    res = part1 * part2 * -1;
    return res;
  }
}
double Package::Check_Anti_T_2(SixJ_coeff &sixJ, SingleP &SP, int J2, int T2,
                               int flag, int index_f, int index_i) {
  double res = 0.0;
  int ch = TP.FindC(J2, T2);
  int a = TP.state[ch][index_i].a;
  int b = TP.state[ch][index_i].b;
  int c = TP.state[ch][index_i].c;
  int E = TP.state[ch][index_i].E;
  int T2z = TP.state[ch][index_i].T_2z;
  int J2ab = TP.state[ch][index_i].Jab_2;
  int T2ab = TP.state[ch][index_i].Tab_2;

  if (flag == 1) {
    // a <-> b
    int phase_t = SP.state[a].j_2 + SP.state[b].j_2 + J2ab + 2 + T2ab;
    int phase = 1 - 2 * ((phase_t / 2) % 2);
    res = vec[ch][index_f][index_i] * phase * -1;
    return res;
  } else {
    // b <-> c
    int size = TP.state[ch].size();
    double part1 = 0.0;
    int phase_t = SP.state[a].j_2 + SP.state[b].j_2 + SP.state[c].j_2 + 3;
    int phase = 1 - 2 * ((phase_t) % 2);
    double part1_t = (J2ab + 1) * (T2ab + 1);
    part1 = phase * std::sqrt(part1_t);
    double part2 = 0.0;
    for (int i = 0; i < size; i++) {
      double part2_t = 0.0;
      if (TP.state[ch][i].E != E || TP.state[ch][i].T_2z != T2z)
        continue;
      // if(TP.state[ch][i].E != E )continue;
      if (TP.state[ch][i].a != c)
        continue;
      if (TP.state[ch][i].b != b || TP.state[ch][i].c != a)
        continue;
      int J2bc = TP.state[ch][i].Jab_2;
      int T2bc = TP.state[ch][i].Tab_2;
      double part2_tt = (J2bc + 1) * (T2bc + 1);
      // phase_t= J2ac + T2ac;
      // phase =1-2*((phase_t/2)%2);
      part2_t = std::sqrt(part2_tt);
      part2_t = part2_t * sixJ.quick(SP.state[b].j_2, SP.state[a].j_2, J2ab, J2,
                                     SP.state[c].j_2, J2bc);
      part2_t = part2_t * sixJ.quick(1, 1, T2ab, T2, 1, T2bc);
      part2 += part2_t * vec[ch][index_f][i];
      // cout<<"\t"<<index_f<<"\t"<<i<<"\t"<<vec[ch][index_f][i]<<endl;
    }
    res = part1 * part2 * -1;
    return res;
  }
}

double Package::Check_Anti_Tz(SixJ_coeff &sixJ, SingleP &SP, int J2, int T_2z,
                              int flag, int index_f, int index_i) {
  double res = 0.0;
  int ch = TP.FindC_Tz(J2, T_2z);
  int a = TP.state_Tz[ch][index_i].a;
  int b = TP.state_Tz[ch][index_i].b;
  int c = TP.state_Tz[ch][index_i].c;
  int E = TP.state_Tz[ch][index_i].E;
  int T2z = TP.state_Tz[ch][index_i].T_2z;
  int J2ab = TP.state_Tz[ch][index_i].Jab_2;
  int T2ab_z = TP.state_Tz[ch][index_i].Tab_2z;

  if (flag == 1) {
    // a <-> b
    int phase_t = SP.state_Tz[a].j_2 + SP.state_Tz[b].j_2 + J2ab;
    int phase = 1 - 2 * ((phase_t / 2) % 2);
    res = vec_Tz[ch][index_f][index_i] * phase;
    if (T2ab_z == 0)
      res = -1 * res;
    return res;
  } else {
    // b <-> c
    int size = TP.state_Tz[ch].size();
    double part1 = 0.0;
    int phase_t = SP.state_Tz[b].j_2 + SP.state_Tz[c].j_2 + J2ab;
    int phase = 1 - 2 * ((phase_t / 2) % 2);
    double part1_t = (J2ab + 1);
    part1 = phase * std::sqrt(part1_t);
    double part2 = 0.0;
    for (int i = 0; i < size; i++) {
      double part2_t = 0.0;
      if (TP.state_Tz[ch][i].E != E || TP.state_Tz[ch][i].T_2z != T2z)
        continue;
      // if(TP.state_Tz[ch][i].E != E )continue;
      if (TP.state_Tz[ch][i].a != a)
        continue;
      if (TP.state_Tz[ch][i].b != c || TP.state_Tz[ch][i].c != b)
        continue;
      int J2ac = TP.state_Tz[ch][i].Jab_2;
      double part2_tt = (J2ac + 1);
      phase_t = J2ac;
      phase = 1 - 2 * ((phase_t / 2) % 2);
      part2_t = phase * std::sqrt(part2_tt);
      part2_t = part2_t * sixJ.quick(SP.state_Tz[a].j_2, SP.state_Tz[b].j_2,
                                     J2ab, J2, SP.state_Tz[c].j_2, J2ac);
      part2 += part2_t * vec_Tz[ch][index_f][i];
    }

    res = part1 * part2;
    if (SP.state_Tz[b].tz_2 * SP.state_Tz[b].tz_2 > 0)
      res = res * (-1); // pp nn

    return res;
  }
}

int Package::FindTPx(int ch, int T2ab, int J2ab, int T2z, int a, int b, int c) {
  int size = TP.state[ch].size();
  int res = -1;
  for (int i = 0; i < size; i++) {
    if (T2ab != TP.state[ch][i].Tab_2 || J2ab != TP.state[ch][i].Jab_2)
      continue;
    if (T2z != TP.state[ch][i].T_2z || a != TP.state[ch][i].a)
      continue;
    if (b != TP.state[ch][i].b || c != TP.state[ch][i].c)
      continue;
    res = i;
    break;
  }
  return res;
  /*
  if(res<0){
    cout<<"channel : "<<ch<<"\t T2ab "<<T2ab<<"\t J2ab "<<J2ab\
        <<"\t T2z "<<T2z<<"\t a "<<a<<"\t b "<<b<<"\t c "<<c<<endl;
    cout<<"wrong at Package::FindTPx"<<endl;
    return res;
    //exit(0);
  }else{
    return res;
  }
  */
}

int Package::FindTPx_Tz(int J2, int T_2z, int J2ab, int a, int b, int c) {
  int ch = TP.FindC_Tz(J2, T_2z);
  int size = TP.state_Tz[ch].size();
  int res = -1;
  for (int i = 0; i < size; i++) {
    if (J2ab != TP.state_Tz[ch][i].Jab_2)
      continue;
    if (a != TP.state_Tz[ch][i].a)
      continue;
    if (b != TP.state_Tz[ch][i].b || c != TP.state_Tz[ch][i].c)
      continue;
    res = i;
    break;
  }
  return res;
  /*
  if(res<0){
    cout<<"channel : "<<ch<<"\t T2ab "<<T2ab<<"\t J2ab "<<J2ab\
        <<"\t T2z "<<T2z<<"\t a "<<a<<"\t b "<<b<<"\t c "<<c<<endl;
    cout<<"wrong at Package::FindTPx"<<endl;
    return res;
    //exit(0);
  }else{
    return res;
  }
  */
}

void Package::build_Tz(CG_coeff &cg, SingleP &SP) {
  int ch_size = TP.channel_Tz.size();
  vec_Tz.resize(ch_size);
  vec_flag_Tz.resize(ch_size);
  // ch_size = 1;
  for (int ch = 0; ch < ch_size; ch++) {
    // cout << "\t --- ch_Tz : " << ch << " --- " << endl;
    int size = TP.state_Tz[ch].size();
    // size = 3;
    // cout << " == ch_size_TP_Tz = " << ch << "\t size = " << size << endl;
    vec_Tz[ch].resize(size);
    vec_flag_Tz[ch].resize(size);
    for (int i = 0; i < size; i++) {
      vec_Tz[ch][i].resize(size);
      vec_flag_Tz[ch][i].resize(size);
      for (int j = 0; j < size; j++) {
        vec_Tz[ch][i][j] = 0.0;
        vec_flag_Tz[ch][i][j] = 0;
      }
    }
    int J2 = TP.channel_Tz[ch].J_2;
    for (int i = 0; i < size; i++) {
      int F_flag = TP.state_Tz[ch][i].F;
      if (F_flag < 0)
        continue;
      int T2ab_z = TP.state_Tz[ch][i].Tab_2z;
      int T2z_i = TP.state_Tz[ch][i].T_2z;
      int J2ab = TP.state_Tz[ch][i].Jab_2;
      // int J2abc = TP.state_Tz[ch][i].J_2;

      int a = TP.state_Tz[ch][i].a;
      int b = TP.state_Tz[ch][i].b;
      int c = TP.state_Tz[ch][i].c;
      int t2z_a = SP.state_Tz[a].tz_2;
      int t2z_b = SP.state_Tz[b].tz_2;
      int t2z_c = SP.state_Tz[c].tz_2;
      int a_x = std::floor(a / 2.0);
      int b_x = std::floor(b / 2.0);
      int c_x = std::floor(c / 2.0);
      for (int j = 0; j < size; j++) {
        int I_flag = TP.state_Tz[ch][j].I;
        if (I_flag < 0)
          continue;

        int Flag = this->TP.FI_compare(F_flag, I_flag);

        if (Flag < 0)
          continue;

        int T2de_z = TP.state_Tz[ch][j].Tab_2z;
        int T2z_j = TP.state_Tz[ch][j].T_2z;
        int J2de = TP.state_Tz[ch][j].Jab_2;
        // int J2def = TP.state_Tz[ch][j].J_2;

        int d = TP.state_Tz[ch][j].a;
        int e = TP.state_Tz[ch][j].b;
        int f = TP.state_Tz[ch][j].c;
        int t2z_d = SP.state_Tz[d].tz_2;
        int t2z_e = SP.state_Tz[e].tz_2;
        int t2z_f = SP.state_Tz[f].tz_2;
        int d_y = std::floor(d / 2.0);
        int e_y = std::floor(e / 2.0);
        int f_y = std::floor(f / 2.0);
        /*cout << ch << "\t" << i << "\t" << j << endl;
        cout << "\t a = " << a << "\t a_x = " << a_x << "\t\t b = " << b
             << "\t b_x = " << b_x << "\t\t c = " << c << "\t c_x = " << c_x
             << endl;
        cout << "\t T2ab_z = " << T2ab_z << "\t t2z_c = " << t2z_c << endl;

        cout << "\t d = " << d << "\t d_y = " << d_y << "\t\t e = " << e
             << "\t e_y = " << e_y << "\t\t f = " << f << "\t f_y = " << f_y
             << endl;
        cout << "\t T2de_z = " << T2de_z << "\t t2z_f = " << t2z_f << endl;
        */
        if (T2z_i != T2z_j)
          continue;
        int T2z = T2z_i;

        double mat = 0.0;
        for (int T2 = 1; T2 <= 3; T2 = T2 + 2) {
          int ch_T = TP.FindC(J2, T2);
          // cout << "T2 = " << T2 << "\t T2z = " << T2z << endl;
          if (std::abs(T2z) > T2)
            continue;
          for (int T2ab = 0; T2ab <= 2; T2ab = T2ab + 2) {
            if (std::abs(T2ab_z) > T2ab)
              continue;
            int x = FindTPx(ch_T, T2ab, J2ab, 0, a_x, b_x, c_x);
            if (x < 0)
              continue;
            double W_x = cg.quick(T2ab, 1, T2, T2ab_z, t2z_c, T2z);
            W_x = W_x * cg.quick(1, 1, T2ab, t2z_a, t2z_b, T2ab_z);
            if (W_x == 0.0)
              continue;
            for (int T2de = 0; T2de <= 2; T2de = T2de + 2) {
              if (std::abs(T2de_z) > T2de)
                continue;
              int y = FindTPx(ch_T, T2de, J2de, 0, d_y, e_y, f_y);
              if (y < 0)
                continue;
              double W_y = cg.quick(T2de, 1, T2, T2de_z, t2z_f, T2z);
              W_y = W_y * cg.quick(1, 1, T2de, t2z_d, t2z_e, T2de_z);
              if (W_y == 0.0)
                continue;
              if (this->vec_flag[ch_T][x][y] != 1) {
                cout << "Wrong happened @ Package::build_Tz" << endl;
                cout << "this->vec_flag[ch_T][x][y] != 1" << endl;
                cout << "ch_T : " << ch_T << "\t x : " << x << "\t y : " << y
                     << endl;
                cout << "ch : " << ch << "\t i : " << i << "\t j : " << j
                     << endl;
                exit(0);
              }
              mat += W_x * W_y * vec[ch_T][x][y];
              /*cout << "\t ch_T = " << ch_T << "\t T_2 = " << T2
                   << "\t T2ab = " << T2ab << "\t T2de = " << T2de
                   << "\t x = " << x << "\t y = " << y
                   << "\t vec = " << vec[ch_T][x][y] << endl;
              cout << "\t W_x = " << W_x << "\t W_y = " << W_y << endl;
             */
            }
          }
        }
        // cout << "mat = " << mat << endl;
        vec_Tz[ch][i][j] = mat;
        vec_flag_Tz[ch][i][j] = 1;
      }
    }
  }
}

void Package::build_S_Tz(CG_coeff &cg, SingleP &SP, SMatrix_Ch &vec_in,
                         SMatrix_Ch &vec_flag_in, SMatrix_Ch &vec_out,
                         SMatrix_Ch &vec_flag_out) {
  cout << "\t Start to build_S_Tz " << endl;
  vec_out.clear();
  vec_flag_out.clear();
  int ch_size = TP.channel_Tz.size();
  vec_out.cal_init(ch_size, 1);
  vec_flag_out.cal_init(ch_size, 0);
  // ch_size = 1;
  for (int ch = 0; ch < ch_size; ch++) {
    cout << "\t --- ch_Tz : " << ch << "\t TP.TPconf_Tz_S.tot_num_ch[ch] : "
         << TP.TPconf_Tz_S.tot_num_ch[ch] << endl;
    int size = TP.state_Tz[ch].size();
    vec_out.cal_beg(size, size);
    vec_flag_out.cal_beg(size, size);

    int J2 = TP.channel_Tz[ch].J_2;
    // int size_conf = TP.TPconf_Tz[ch].size();
    // vector<vector<double>> vec_t;
    // vector<vector<int>> vec_flag_t;
    // vec_t.resize(size);
    // vec_flag_t.resize(size);
    // for (int i = 0; i < size; i++) {
    //   vec_t[i].resize(size);
    //   vec_flag_t[i].resize(size);
    //   for (int j = 0; j < size; j++) {
    //     vec_t[i][j] = 0.0;
    //     vec_flag_t[i][j] = 0;
    //   }
    // }
    // for (int co = 0; co < size_conf; co++) {
    //   int i = TP.TPconf_Tz[ch][co].F;
    //   int j = TP.TPconf_Tz[ch][co].I;
    for (int i = 0; i < size; i++) {
      vec_out.cal_row_off(ch, i);
      vec_flag_out.cal_row_off(ch, i);
      int F_flag = TP.state_Tz[ch][i].F;
      if (F_flag < 0)
        continue;
      int T2ab_z = TP.state_Tz[ch][i].Tab_2z;
      int T2z_i = TP.state_Tz[ch][i].T_2z;
      int J2ab = TP.state_Tz[ch][i].Jab_2;
      // int J2abc = TP.state_Tz[ch][i].J_2;

      int a = TP.state_Tz[ch][i].a;
      int b = TP.state_Tz[ch][i].b;
      int c = TP.state_Tz[ch][i].c;
      int t2z_a = SP.state_Tz[a].tz_2;
      int t2z_b = SP.state_Tz[b].tz_2;
      int t2z_c = SP.state_Tz[c].tz_2;
      int a_x = std::floor(a / 2.0);
      int b_x = std::floor(b / 2.0);
      int c_x = std::floor(c / 2.0);

      for (int j = 0; j < size; j++) {
        int I_flag = TP.state_Tz[ch][j].I;
        if (I_flag < 0)
          continue;
        int Flag = this->TP.FI_compare(F_flag, I_flag);
        if (Flag < 0)
          continue;
        if (TP.TPconf_Tz_S.val(ch, i, j) != 1)
          continue;

        int T2de_z = TP.state_Tz[ch][j].Tab_2z;
        int T2z_j = TP.state_Tz[ch][j].T_2z;
        int J2de = TP.state_Tz[ch][j].Jab_2;
        // int J2def = TP.state_Tz[ch][j].J_2;

        int d = TP.state_Tz[ch][j].a;
        int e = TP.state_Tz[ch][j].b;
        int f = TP.state_Tz[ch][j].c;
        int t2z_d = SP.state_Tz[d].tz_2;
        int t2z_e = SP.state_Tz[e].tz_2;
        int t2z_f = SP.state_Tz[f].tz_2;
        int d_y = std::floor(d / 2.0);
        int e_y = std::floor(e / 2.0);
        int f_y = std::floor(f / 2.0);
        /*cout << ch << "\t" << i << "\t" << j << endl;
        cout << "\t a = " << a << "\t a_x = " << a_x << "\t\t b = " << b
             << "\t b_x = " << b_x << "\t\t c = " << c << "\t c_x = " << c_x
             << endl;
        cout << "\t T2ab_z = " << T2ab_z << "\t t2z_c = " << t2z_c << endl;

        cout << "\t d = " << d << "\t d_y = " << d_y << "\t\t e = " << e
             << "\t e_y = " << e_y << "\t\t f = " << f << "\t f_y = " << f_y
             << endl;
        cout << "\t T2de_z = " << T2de_z << "\t t2z_f = " << t2z_f << endl;
        */
        if (T2z_i != T2z_j)
          continue;
        int T2z = T2z_i;

        double mat = 0.0;
        for (int T2 = 1; T2 <= 3; T2 = T2 + 2) {
          int ch_T = TP.FindC(J2, T2);
          // cout << "T2 = " << T2 << "\t T2z = " << T2z << endl;
          if (std::abs(T2z) > T2)
            continue;
          for (int T2ab = 0; T2ab <= 2; T2ab = T2ab + 2) {
            if (std::abs(T2ab_z) > T2ab)
              continue;
            int x = FindTPx(ch_T, T2ab, J2ab, 0, a_x, b_x, c_x);
            if (x < 0)
              continue;
            double W_x = cg.quick(T2ab, 1, T2, T2ab_z, t2z_c, T2z);
            W_x = W_x * cg.quick(1, 1, T2ab, t2z_a, t2z_b, T2ab_z);
            if (W_x == 0.0)
              continue;
            for (int T2de = 0; T2de <= 2; T2de = T2de + 2) {
              if (std::abs(T2de_z) > T2de)
                continue;
              int y = FindTPx(ch_T, T2de, J2de, 0, d_y, e_y, f_y);
              if (y < 0)
                continue;
              double W_y = cg.quick(T2de, 1, T2, T2de_z, t2z_f, T2z);
              W_y = W_y * cg.quick(1, 1, T2de, t2z_d, t2z_e, T2de_z);
              if (W_y == 0.0)
                continue;
              // if (this->vec_flag[ch_T][x][y] != 1) {
              if (vec_flag_in.val(ch_T, x, y) != 1) {
                cout << "Wrong happened @ Package::build_Tz" << endl;
                cout << "this->vec_flag[ch_T][x][y] != 1" << endl;
                cout << "ch_T : " << ch_T << "\t x : " << x << "\t y : " << y
                     << endl;
                cout << "ch : " << ch << "\t i : " << i << "\t j : " << j
                     << endl;
                exit(0);
              }
              // mat += W_x * W_y * vec[ch_T][x][y];
              mat += W_x * W_y * vec_in.val(ch_T, x, y);
              /*cout << "\t ch_T = " << ch_T << "\t T_2 = " << T2
                   << "\t T2ab = " << T2ab << "\t T2de = " << T2de
                   << "\t x = " << x << "\t y = " << y
                   << "\t vec = " << vec[ch_T][x][y] << endl;
              cout << "\t W_x = " << W_x << "\t W_y = " << W_y << endl;
             */
            }
          }
        }
        // cout << "mat = " << mat << endl;
        // vec_Tz[ch][i][j] = mat;
        // vec_flag_Tz[ch][i][j] = 1;
        vec_out.cal_col_ind(ch, i, j, mat);
        vec_flag_out.cal_col_ind(ch, i, j, 1);
      }
      // vec_t[i][j] = mat;
      // vec_flag_t[i][j] = 1;
    }
    // SMatrix vec_t_S;
    // SMatrix vec_flag_t_S;
    // vec_t_S.build(vec_t);
    // cout << " -- finished vec_t_S.build(vec_t) -- " << endl;
    // vec_flag_t_S.build(vec_flag_t);
    // cout << " -- vec_flag_t_S.build(vec_flag_t) -- " << endl;
    // vec_out.row_num_ch[ch] = vec_t_S.row_num;
    // vec_out.col_num_ch[ch] = vec_t_S.col_num;
    // vec_out.tot_num_ch[ch] = vec_t_S.tot_num;
    // vec_out.row_off_ch[ch] = vec_t_S.row_offsets;
    // vec_out.col_ind_ch[ch] = vec_t_S.col_indices;
    // vec_out.vec_D_ch[ch] = vec_t_S.vec_D;
    //
    // vec_flag_out.row_num_ch[ch] = vec_flag_t_S.row_num;
    // vec_flag_out.col_num_ch[ch] = vec_flag_t_S.col_num;
    // vec_flag_out.tot_num_ch[ch] = vec_flag_t_S.tot_num;
    // vec_flag_out.row_off_ch[ch] = vec_flag_t_S.row_offsets;
    // vec_flag_out.col_ind_ch[ch] = vec_flag_t_S.col_indices;
    // vec_flag_out.vec_I_ch[ch] = vec_flag_t_S.vec_I;
    //
    // vec_t_S.clear();
    // vec_flag_t_S.clear();
    // vec_t.clear();
    // vec_flag_t.clear();
    vec_out.cal_end(ch);
    vec_flag_out.cal_end(ch);
  }
}

void Package::build_TzM(CG_coeff &cg, SingleP &SP) {
  int size = TP.state_TzM.size();
  vec_TzM.resize(size);
  int J2max = TP.J2_max;

  for (int i = 0; i < size; i++) {
    vec_TzM[i].resize(size);
    for (int j = 0; j < size; j++) {
      vec_TzM[i][j] = 0.0;
    }
  }

  for (int i = 0; i < size; i++) {
    int T2z_i = TP.state_TzM[i].T_2z;
    int J2m_i = TP.state_TzM[i].J_2m;
    int T2ab_z = TP.state_TzM[i].Tab_2z;
    int J2ab_m = TP.state_TzM[i].Jab_2m;
    int a_m = TP.state_TzM[i].a_m;
    int b_m = TP.state_TzM[i].b_m;
    int c_m = TP.state_TzM[i].c_m;
    int j2_a = SP.state_TzM[a_m].j_2;
    int j2_b = SP.state_TzM[b_m].j_2;
    int j2_c = SP.state_TzM[c_m].j_2;
    int m2_a = SP.state_TzM[a_m].m_2;
    int m2_b = SP.state_TzM[b_m].m_2;
    int m2_c = SP.state_TzM[c_m].m_2;
    int a = SP.index_M2J[a_m];
    int b = SP.index_M2J[b_m];
    int c = SP.index_M2J[c_m];

    for (int j = 0; j < size; j++) {
      int T2z_j = TP.state_TzM[j].T_2z;
      int J2m_j = TP.state_TzM[j].J_2m;
      int T2de_z = TP.state_TzM[j].Tab_2z;
      int J2de_m = TP.state_TzM[j].Jab_2m;
      int d_m = TP.state_TzM[j].a_m;
      int e_m = TP.state_TzM[j].b_m;
      int f_m = TP.state_TzM[j].c_m;
      int j2_d = SP.state_TzM[d_m].j_2;
      int j2_e = SP.state_TzM[e_m].j_2;
      int j2_f = SP.state_TzM[f_m].j_2;
      int m2_d = SP.state_TzM[d_m].m_2;
      int m2_e = SP.state_TzM[e_m].m_2;
      int m2_f = SP.state_TzM[f_m].m_2;
      int d = SP.index_M2J[d_m];
      int e = SP.index_M2J[e_m];
      int f = SP.index_M2J[f_m];
      // cout<<i<<"\t"<<j<<endl;
      if (T2z_i != T2z_j)
        continue;
      int T2z = T2z_j;
      double mat = 0.0;
      int J2_min = std::max(std::abs(J2m_i), std::abs(J2m_j));
      for (int J2 = J2_min; J2 <= J2max; J2 = J2 + 2) {
        int ch_Tz = TP.FindC_Tz(J2, T2z);
        for (int J2ab = std::abs(J2ab_m); J2ab <= J2max; J2ab = J2ab + 2) {
          int x = FindTPx_Tz(J2, T2z, J2ab, a, b, c);
          if (x < 0)
            continue;
          double W_x = cg.quick(J2ab, j2_c, J2, J2ab_m, m2_c, J2m_i);
          W_x = W_x * cg.quick(j2_a, j2_b, J2ab, m2_a, m2_b, J2ab_m);
          if (W_x == 0.0)
            continue;

          for (int J2de = std::abs(J2de_m); J2de <= J2max; J2de = J2de + 2) {
            int y = FindTPx_Tz(J2, T2z, J2de, d, e, f);
            if (y < 0)
              continue;
            double W_y = cg.quick(J2de, j2_f, J2, J2de_m, m2_f, J2m_j);
            W_y = W_y * cg.quick(j2_d, j2_e, J2de, m2_d, m2_e, J2de_m);
            if (W_y == 0.0)
              continue;
            mat += W_x * W_y * vec_Tz[ch_Tz][x][y];
          }
        }
      }
      vec_TzM[i][j] = mat;
    }
  }
}

void Package::build_anti_bc(SixJ_coeff &sixJ, SingleP &SP) {
  int ch_size = vec.size();
  vec_anti.resize(ch_size);
  for (int ch = 0; ch < ch_size; ch++) {
    int size = vec[ch].size();
    int J2 = TP.channel[ch].J_2;
    int T2 = TP.channel[ch].T_2;
    vec_anti[ch].resize(size);
    // cout<<"Channel : "<<ch<<endl;
    for (int i = 0; i < size; i++) {
      vec_anti[ch][i].resize(size);
      for (int j = 0; j < size; j++) {
        double res = Check_Anti_T_2(sixJ, SP, J2, T2, 2, i, j);
        vec_anti[ch][i][j] = res;
        if (std::abs(res) > 0.000001) {
          // cout<<i<<"\t"<<j<<"\t"<<res<<endl;
        }
      }
    }
  }
}

void Package::build_anti_ab_Tz(SixJ_coeff &sixJ, SingleP &SP) {
  int ch_size = vec_Tz.size();
  vec_anti_Tz.resize(ch_size);
  for (int ch = 0; ch < ch_size; ch++) {
    int size = vec_Tz[ch].size();
    int J2 = TP.channel_Tz[ch].J_2;
    int T_2z = TP.channel_Tz[ch].T_2z;
    vec_anti_Tz[ch].resize(size);
    // cout<<"Channel : "<<ch<<endl;
    for (int i = 0; i < size; i++) {
      vec_anti_Tz[ch][i].resize(size);
      for (int j = 0; j < size; j++) {
        double res = Check_Anti_Tz(sixJ, SP, J2, T_2z, 1, i, j);
        vec_anti_Tz[ch][i][j] = res;
        if (std::abs(res) > 0.000001) {
          // cout<<i<<"\t"<<j<<"\t"<<res<<endl;
        }
      }
    }
  }
}
void Package::build_anti_bc_Tz(SixJ_coeff &sixJ, SingleP &SP) {
  int ch_size = vec_Tz.size();
  vec_anti_Tz.resize(ch_size);
  for (int ch = 0; ch < ch_size; ch++) {
    int size = vec_Tz[ch].size();
    int J2 = TP.channel_Tz[ch].J_2;
    int T_2z = TP.channel_Tz[ch].T_2z;
    vec_anti_Tz[ch].resize(size);
    // cout<<"Channel : "<<ch<<endl;
    for (int i = 0; i < size; i++) {
      vec_anti_Tz[ch][i].resize(size);
      for (int j = 0; j < size; j++) {
        double res = Check_Anti_Tz(sixJ, SP, J2, T_2z, 2, i, j);
        vec_anti_Tz[ch][i][j] = res;
        if (std::abs(res) > 0.000001) {
          // cout<<i<<"\t"<<j<<"\t"<<res<<endl;
        }
      }
    }
  }
}

void Package::print_nnlosat(SingleP SP_t, std::string sp_nnlosat_file,
                            std::string filename) {
  // vector<State_SP> sp_nnlosat_state;
  // vector<int> z2sat;
  sp_nnlosat_state.clear();
  z2sat.clear();

  std::ifstream sp_file;
  sp_file.open(sp_nnlosat_file);
  char u1[256]; // unuse
  sp_file.getline(u1, 256, '\n');
  while (!sp_file.eof()) {
    State_SP sp_nnlosat_t;
    int i, N, n, l, j2;
    sp_file >> i >> N >> n >> l >> j2;
    sp_nnlosat_t.index = i;
    sp_nnlosat_t.E = N;
    sp_nnlosat_t.n = n;
    sp_nnlosat_t.l = l;
    sp_nnlosat_t.j_2 = j2;
    sp_nnlosat_state.push_back(sp_nnlosat_t);
  }
  int size_sp = SP_t.state.size();
  z2sat.resize(size_sp);
  for (int i = 0; i < size_sp; i++) {
    int n = SP_t.state[i].n;
    int l = SP_t.state[i].l;
    int j2 = SP_t.state[i].j_2;
    int E = SP_t.state[i].E;
    int flag = -1;
    z2sat[i] = -1;
    for (int j = 0; j < sp_nnlosat_state.size(); j++) {
      int n_sat = sp_nnlosat_state[j].n;
      int l_sat = sp_nnlosat_state[j].l;
      int E_sat = sp_nnlosat_state[j].E;
      int j2_sat = sp_nnlosat_state[j].j_2;
      if (n != n_sat || l != l_sat)
        continue;
      if (j2 != j2_sat || E != E_sat)
        continue;
      z2sat[i] = sp_nnlosat_state[j].index;
      flag = 1;
    }
    if (flag == -1) {
      cout << "Wrong happend at Package::print_nnlosat " << endl;
      cout << "Can't find SP from SP_nnlosat : i = " << i << endl;
      // exit(0);
    }
  }
  for (int i = 0; i < z2sat.size(); i++) {
    cout << "sp : " << i << "\t sp_nnlosat i :" << z2sat[i] << endl;
  }

  std::ofstream file;
  file.open(filename);
  int size_ch = vec.size();
  file
      << "i\t j\t k\t l\t m\t n\t 2ij\t 2kl\t 2J\t Tij\t Tlm\t 2T\t <ijk|V|lmn>"
      << endl;
  for (int ch = 0; ch < size_ch; ch++) {
    int size = vec[ch].size();
    // file << "Channel : " << ch << endl;
    for (int i = 0; i < size; i++) {
      int a_p = TP.state[ch][i].a;
      int b_p = TP.state[ch][i].b;
      int c_p = TP.state[ch][i].c;
      int J2ab_p = TP.state[ch][i].Jab_2;
      int Tab_p = TP.state[ch][i].Tab_2 / 2;
      int J2_p = TP.state[ch][i].J_2;
      int T2_p = TP.state[ch][i].T_2;
      int ix = z2sat[a_p];
      int jx = z2sat[b_p];
      int kx = z2sat[c_p];

      for (int j = 0; j < size; j++) {
        int a = TP.state[ch][j].a;
        int b = TP.state[ch][j].b;
        int c = TP.state[ch][j].c;
        int J2ab = TP.state[ch][j].Jab_2;
        int Tab = TP.state[ch][j].Tab_2 / 2;
        int J2 = TP.state[ch][j].J_2;
        int T2 = TP.state[ch][j].T_2;

        int lx = z2sat[a];
        int mx = z2sat[b];
        int nx = z2sat[c];

        if (std::abs(vec[ch][i][j]) > 0.0001) {
          file << ix << "\t" << jx << "\t" << kx << "\t" << lx << "\t" << mx
               << "\t" << nx << "\t" << J2ab_p << "\t" << J2ab << "\t" << J2
               << "\t" << Tab_p << "\t" << Tab << "\t" << T2 << "\t\t"
               << vec[ch][i][j] << "\t" << i << "\t" << j << endl;
        }
      }
    }
  }
}

void Package::Readnnlosat(SingleP SP_t, std::string sp_nnlosat_file,
                          std::string me_nnlosat_file, vector<int> abc_range) {

  int a_min = abc_range[0];
  int a_max = abc_range[1];
  int b_min = abc_range[2];
  int b_max = abc_range[3];
  int c_min = abc_range[4];
  int c_max = abc_range[5];

  sp_nnlosat_state.clear();
  z2sat.clear();
  sat2z.clear();
  // vector<State_SP> sp_nnlosat_state;
  // vector<int> z2sat;
  // vector<int> sat2z;

  std::ifstream sp_file;
  sp_file.open(sp_nnlosat_file);
  char u1[256]; // unuse
  sp_file.getline(u1, 256, '\n');
  while (!sp_file.eof()) {
    State_SP sp_nnlosat_t;
    int i, N, n, l, j2;
    sp_file >> i >> N >> n >> l >> j2;
    sp_nnlosat_t.index = i;
    sp_nnlosat_t.E = N;
    sp_nnlosat_t.n = n;
    sp_nnlosat_t.l = l;
    sp_nnlosat_t.j_2 = j2;
    sp_nnlosat_state.push_back(sp_nnlosat_t);
  }
  int size_sp = SP_t.state.size();
  z2sat.resize(size_sp);
  for (int i = 0; i < size_sp; i++) {
    int n = SP_t.state[i].n;
    int l = SP_t.state[i].l;
    int j2 = SP_t.state[i].j_2;
    int E = SP_t.state[i].E;
    int flag = -1;
    z2sat[i] = -1;
    for (int j = 0; j < sp_nnlosat_state.size(); j++) {
      int n_sat = sp_nnlosat_state[j].n;
      int l_sat = sp_nnlosat_state[j].l;
      int E_sat = sp_nnlosat_state[j].E;
      int j2_sat = sp_nnlosat_state[j].j_2;
      if (n != n_sat || l != l_sat)
        continue;
      if (j2 != j2_sat || E != E_sat)
        continue;
      z2sat[i] = sp_nnlosat_state[j].index;
      flag = 1;
    }
    if (flag == -1) {
      cout << "Wrong happend at Package::print_nnlosat " << endl;
      cout << "Can't find SP from SP_nnlosat : i = " << i << endl;
      // exit(0);
    }
  }
  for (int i = 0; i < z2sat.size(); i++) {
    // cout << "sp : " << i << "\t sp_nnlosat i :" << z2sat[i] << endl;
  }
  sat2z.resize(sp_nnlosat_state.size());
  for (int j = 0; j < sp_nnlosat_state.size(); j++) {
    int n_sat = sp_nnlosat_state[j].n;
    int l_sat = sp_nnlosat_state[j].l;
    int E_sat = sp_nnlosat_state[j].E;
    int j2_sat = sp_nnlosat_state[j].j_2;
    int flag = -1;
    sat2z[j] = -1;
    for (int i = 0; i < size_sp; i++) {
      int n = SP_t.state[i].n;
      int l = SP_t.state[i].l;
      int j2 = SP_t.state[i].j_2;
      int E = SP_t.state[i].E;

      if (n != n_sat || l != l_sat)
        continue;
      if (j2 != j2_sat || E != E_sat)
        continue;
      sat2z[j] = i;
      flag = 1;
    }
    if (flag == -1) {
      cout << "Wrong happend at compareVSnnlosat " << endl;
      cout << "Can't find SP_nnlosat from SP : j = " << j << endl;
      // exit(0);
    }
  }
  for (int j = 0; j < sat2z.size(); j++) {
    cout << "sp_nnlosat : " << j << "\t j+1 " << j + 1
         << "\t sp i :" << sat2z[j] << endl;
  }

  int size_ch_TP = this->vec.size();

  this->nnlosat_vec.resize(size_ch_TP);
  for (int ch = 0; ch < size_ch_TP; ch++) {
    int size_i = this->vec[ch].size();
    nnlosat_vec[ch].resize(size_i);
    for (int j = 0; j < size_i; j++) {
      nnlosat_vec[ch][j].resize(size_i);
    }
  }
  std::ifstream me_file;
  me_file.open(me_nnlosat_file);
  me_file.getline(u1, 256, '\n');
  me_file.getline(u1, 256, '\n');
  while (!me_file.eof()) {
    int i, j, k, l, m, n, J2ij, J2lm, J2, Tij, Tlm, T2;
    double val;
    me_file >> i >> j >> k >> l >> m >> n >> J2ij >> J2lm >> J2 >> Tij >> Tlm >>
        T2 >> val;
    int sp_size = sat2z.size();
    if (i > sp_size || j > sp_size || k > sp_size)
      continue;
    if (l > sp_size || m > sp_size || n > sp_size)
      continue;
    int a_p = sat2z[i - 1];
    int b_p = sat2z[j - 1];
    int c_p = sat2z[k - 1];
    int J2ab_p = J2ij;
    int T2ab_p = Tij * 2;

    int a = sat2z[l - 1];
    int b = sat2z[m - 1];
    int c = sat2z[n - 1];
    int J2ab = J2lm;
    int T2ab = Tlm * 2;

    if (a_p < 0 || b_p < 0 || c_p < 0)
      continue;
    if (a < 0 || b < 0 || c < 0)
      continue;

    if (a_p < a_min || a_p > a_max)
      continue;
    if (b_p < b_min || b_p > b_max)
      continue;
    if (c_p < c_min || c_p > c_max)
      continue;

    if (a < a_min || a > a_max)
      continue;
    if (b < b_min || b > b_max)
      continue;
    if (c < c_min || c > c_max)
      continue;

    int ch_TP = TP.FindC(J2, T2);
    int tp_final = FindTP_index(a_p, b_p, c_p, J2ab_p, T2ab_p, J2, T2);
    int tp_init = FindTP_index(a, b, c, J2ab, T2ab, J2, T2);
    if (tp_final < 0 || tp_init < 0) {
      cout
          << "Wrong happend at Package::Readnnlosat tp_final < 0 or tp_init < 0"
          << endl;
      cout << "tp_final = " << tp_final << "\t tp_init = " << tp_init << endl;

      cout << "i\t j\t k\t l\t m\t n\t 2ij\t 2kl\t 2J\t Tij\t Tlm\t 2T\t "
              "<ijk|V|lmn>"
           << endl;
      cout << i << "\t" << j << "\t" << k << "\t" << l << "\t" << m << "\t" << n
           << "\t" << J2ij << "\t" << J2lm << "\t" << J2 << "\t" << Tij << "\t"
           << Tlm << "\t" << T2 << "\t" << val << endl;

      cout << "a_p\t b_p\t c_p\t a\t b\t c\t J2ab_p\t J2ab\t T2ab_p\t T2ab\t "
              "J2\t T2\t "
           << endl;
      cout << a_p << "\t" << b_p << "\t" << c_p << "\t" << a << "\t" << b
           << "\t" << c << "\t" << J2ab_p << "\t" << J2ab << "\t" << T2ab_p
           << "\t" << T2ab << "\t" << J2 << "\t" << T2 << "\t " << val << endl;

      exit(0);
      continue;
    }
    nnlosat_vec[ch_TP][tp_final][tp_init] = val;
    nnlosat_vec[ch_TP][tp_init][tp_final] = val;
  }
  me_file.close();
}

void Package::Readnnlosat_SP(SingleP SP_t, std::string sp_nnlosat_file) {

  sp_nnlosat_state.clear();
  z2sat.clear();
  sat2z.clear();
  // vector<State_SP> sp_nnlosat_state;
  // vector<int> z2sat;
  // vector<int> sat2z;

  std::ifstream sp_file;
  sp_file.open(sp_nnlosat_file);
  char u1[256]; // unuse
  sp_file.getline(u1, 256, '\n');
  while (!sp_file.eof()) {
    State_SP sp_nnlosat_t;
    int i, N, n, l, j2;
    sp_file >> i >> N >> n >> l >> j2;
    sp_nnlosat_t.index = i;
    sp_nnlosat_t.E = N;
    sp_nnlosat_t.n = n;
    sp_nnlosat_t.l = l;
    sp_nnlosat_t.j_2 = j2;
    sp_nnlosat_state.push_back(sp_nnlosat_t);
  }
  int size_sp = SP_t.state.size();
  z2sat.resize(size_sp);
  for (int i = 0; i < size_sp; i++) {
    int n = SP_t.state[i].n;
    int l = SP_t.state[i].l;
    int j2 = SP_t.state[i].j_2;
    int E = SP_t.state[i].E;
    int flag = -1;
    z2sat[i] = -1;
    for (int j = 0; j < sp_nnlosat_state.size(); j++) {
      int n_sat = sp_nnlosat_state[j].n;
      int l_sat = sp_nnlosat_state[j].l;
      int E_sat = sp_nnlosat_state[j].E;
      int j2_sat = sp_nnlosat_state[j].j_2;
      if (n != n_sat || l != l_sat)
        continue;
      if (j2 != j2_sat || E != E_sat)
        continue;
      z2sat[i] = sp_nnlosat_state[j].index;
      flag = 1;
    }
    if (flag == -1) {
      cout << "Wrong happend at Package::print_nnlosat " << endl;
      cout << "Can't find SP from SP_nnlosat : i = " << i << endl;
      // exit(0);
    }
  }
  for (int i = 0; i < z2sat.size(); i++) {
    // cout << "sp : " << i << "\t sp_nnlosat i :" << z2sat[i] << endl;
  }
  sat2z.resize(sp_nnlosat_state.size());
  for (int j = 0; j < sp_nnlosat_state.size(); j++) {
    int n_sat = sp_nnlosat_state[j].n;
    int l_sat = sp_nnlosat_state[j].l;
    int E_sat = sp_nnlosat_state[j].E;
    int j2_sat = sp_nnlosat_state[j].j_2;
    int flag = -1;
    sat2z[j] = -1;
    for (int i = 0; i < size_sp; i++) {
      int n = SP_t.state[i].n;
      int l = SP_t.state[i].l;
      int j2 = SP_t.state[i].j_2;
      int E = SP_t.state[i].E;

      if (n != n_sat || l != l_sat)
        continue;
      if (j2 != j2_sat || E != E_sat)
        continue;
      sat2z[j] = i;
      flag = 1;
    }
    if (flag == -1) {
      cout << "Wrong happend at compareVSnnlosat " << endl;
      cout << "Can't find SP_nnlosat from SP : j = " << j << endl;
      // exit(0);
    }
  }
  for (int j = 0; j < sat2z.size(); j++) {
    cout << "sp_nnlosat : " << j << "\t j+1 " << j + 1
         << "\t sp i :" << sat2z[j] << endl;
  }
}

void Package::CompareVsNNlosat(string filename, string filename2,
                               string filename3) {
  std::ofstream file;
  file.open(filename);
  std::ofstream file2;
  file2.open(filename2);
  std::ofstream file3;
  file3.open(filename3);
  int size_ch = vec.size();
  for (int ch = 0; ch < size_ch; ch++) {
    int size = vec[ch].size();
    file << "Channel : " << ch << endl;
    file << "a' b' c'\ta  b  c  J2ab'  J2ab  J2  T2ab'  T2ab  T2 ";
    file << "\t i\t j\t val_sat\t val_YZ\t\t Diff\t\t ratio" << endl;
    file2 << "Channel : " << ch << endl;
    file2 << "a' b' c'\ta  b  c  J2ab'  J2ab  J2  T2ab'  T2ab  T2 ";
    file2 << "\t i\t j\t val_sat\t val_YZ\t\t Diff\t\t ratio" << endl;
    file3 << "Channel : " << ch << endl;
    file3 << "a' b' c'\ta  b  c  J2ab'  J2ab  J2  T2ab'  T2ab  T2 ";
    file3 << "\t i\t j\t val_sat\t val_YZ\t\t Diff\t\t ratio" << endl;
    for (int i = 0; i < size; i++) {
      for (int j = 0; j <= i; j++) {

        if (std::abs(vec[ch][i][j]) > 0.000001 ||
            std::abs(nnlosat_vec[ch][i][j]) > 0.000001) {
          double diff = nnlosat_vec[ch][i][j] - vec[ch][i][j];
          if (std::abs(diff) < 0.0000001)
            diff = 0.0;
          double ratio = nnlosat_vec[ch][i][j] / vec[ch][i][j];
          int a_p = TP.state[ch][i].a;
          int b_p = TP.state[ch][i].b;
          int c_p = TP.state[ch][i].c;
          int J2ab_p = TP.state[ch][i].Jab_2;
          int T2ab_p = TP.state[ch][i].Tab_2;

          int a = TP.state[ch][j].a;
          int b = TP.state[ch][j].b;
          int c = TP.state[ch][j].c;
          int J2ab = TP.state[ch][j].Jab_2;
          int T2ab = TP.state[ch][j].Tab_2;

          int J2 = TP.state[ch][j].J_2;
          int T2 = TP.state[ch][j].T_2;

          int a_sat = z2sat[a];
          int b_sat = z2sat[b];
          int c_sat = z2sat[c];

          int a_sat_p = z2sat[a_p];
          int b_sat_p = z2sat[b_p];
          int c_sat_p = z2sat[c_p];

          if (a_sat_p < b_sat_p || b_sat_p < c_sat_p || a_sat_p < c_sat_p)
            continue;
          if (a_sat < b_sat || b_sat < c_sat || a_sat < c_sat)
            continue;
          /*
          if (a_p == 2 && b_p == 1 && c_p == 0)
            continue;
          if (a_p == 2 && b_p == 1 && c_p == 1)
            continue;
          if (a_p == 2 && b_p == 2 && c_p == 1)
            continue;

          if (a == 2 && b == 1 && c == 0)
            continue;
          if (a == 2 && b == 1 && c == 1)
            continue;
          if (a == 2 && b == 2 && c == 1)
            continue; */

          file << a_p << "  " << b_p << "  " << c_p << "\t\t" << a << "  " << b
               << "  " << c << "    " << J2ab_p << "     " << J2ab << "    "
               << J2 << "    " << T2ab_p << "      " << T2ab << "    " << T2
               << "\t " << i << "\t" << j << "\t" << std::setw(8)
               << nnlosat_vec[ch][i][j] << "\t " << std::setw(8)
               << vec[ch][i][j] << "\t" << std::setw(8) << diff << "\t"
               << std::setw(8) << ratio;
          if (std::abs(ratio - 1) > 0.02)
            file << "\t * \t " << std::setw(3) << std::abs(ratio - 1) * 100
                 << "%";
          file << endl;

          int flag = 1;
          if (a_sat_p >= b_sat_p && b_sat_p >= c_sat_p) {
            if (a_sat >= b_sat && b_sat >= c_sat) {
              flag = -1;
            }
          }
          if (flag > 0) {
            file2 << a_p << "  " << b_p << "  " << c_p << "\t\t" << a << "  "
                  << b << "  " << c << "    " << J2ab_p << "     " << J2ab
                  << "    " << J2 << "    " << T2ab_p << "      " << T2ab
                  << "    " << T2 << "\t " << i << "\t" << j << "\t"
                  << std::setw(8) << nnlosat_vec[ch][i][j] << "\t "
                  << std::setw(8) << vec[ch][i][j] << "\t" << std::setw(8)
                  << diff << "\t" << std::setw(8) << ratio;
            if (std::abs(ratio - 1) > 0.02)
              file2 << "\t *";
            file2 << endl;
          }
          if ((vec[ch][i][j] * nnlosat_vec[ch][i][j]) < 0) {
            file3 << a_p << "  " << b_p << "  " << c_p << "\t\t" << a << "  "
                  << b << "  " << c << "    " << J2ab_p << "     " << J2ab
                  << "    " << J2 << "    " << T2ab_p << "      " << T2ab
                  << "    " << T2 << "\t " << i << "\t" << j << "\t"
                  << std::setw(8) << nnlosat_vec[ch][i][j] << "\t "
                  << std::setw(8) << vec[ch][i][j] << "\t" << std::setw(8)
                  << diff << "\t" << std::setw(8) << ratio;
            if (std::abs(ratio - 1) > 0.02)
              file3 << "\t *";
            file3 << endl;
          }
        }
      }
    }
  }
  file.close();
  file2.close();
  file3.close();
}

void Package::Compare_satVSzhuo(string me_nnlosat_file, string filename2,
                                vector<int> abc_range) {

  int a_min = abc_range[0];
  int a_max = abc_range[1];
  int b_min = abc_range[2];
  int b_max = abc_range[3];
  int c_min = abc_range[4];
  int c_max = abc_range[5];

  char u1[256]; // unuse

  std::ofstream out_file;
  out_file.open(filename2);
  out_file << "i\t j\t k\t l\t m\t n\t 2ij\t 2kl\t 2J\t Tij\t Tlm\t 2T\t "
              "<ijk|V|lmn>\t";
  out_file << "mat_YZ\t rel diff (output when > 2% )" << endl;
  std::ifstream me_file;
  me_file.open(me_nnlosat_file);
  me_file.getline(u1, 256, '\n');
  me_file.getline(u1, 256, '\n');
  while (!me_file.eof()) {
    int i, j, k, l, m, n, J2ij, J2lm, J2, Tij, Tlm, T2;
    double val;
    me_file >> i >> j >> k >> l >> m >> n >> J2ij >> J2lm >> J2 >> Tij >> Tlm >>
        T2 >> val;
    int sp_size = sat2z.size();
    if (i > sp_size || j > sp_size || k > sp_size)
      continue;
    if (l > sp_size || m > sp_size || n > sp_size)
      continue;
    int a_p = sat2z[i - 1];
    int b_p = sat2z[j - 1];
    int c_p = sat2z[k - 1];
    int J2ab_p = J2ij;
    int T2ab_p = Tij * 2;

    int a = sat2z[l - 1];
    int b = sat2z[m - 1];
    int c = sat2z[n - 1];
    int J2ab = J2lm;
    int T2ab = Tlm * 2;

    if (a_p < 0 || b_p < 0 || c_p < 0)
      continue;
    if (a < 0 || b < 0 || c < 0)
      continue;

    if (a_p < a_min || a_p > a_max)
      continue;
    if (b_p < b_min || b_p > b_max)
      continue;
    if (c_p < c_min || c_p > c_max)
      continue;

    if (a < a_min || a > a_max)
      continue;
    if (b < b_min || b > b_max)
      continue;
    if (c < c_min || c > c_max)
      continue;

    if (J2 > 9) {
      cout << "J2 > 9 continue Becarefule" << endl;
      continue;
    }

    int ch_TP = TP.FindC(J2, T2);
    int tp_final = FindTP_index(a_p, b_p, c_p, J2ab_p, T2ab_p, J2, T2);
    int tp_init = FindTP_index(a, b, c, J2ab, T2ab, J2, T2);
    if (tp_final < 0 || tp_init < 0)
      continue;
    if (tp_final < 0 || tp_init < 0) {
      cout
          << "Wrong happend at Package::Readnnlosat tp_final < 0 or tp_init < 0"
          << endl;
      cout << "tp_final = " << tp_final << "\t tp_init = " << tp_init << endl;

      cout << "i\t j\t k\t l\t m\t n\t 2ij\t 2kl\t 2J\t Tij\t Tlm\t 2T\t "
              "<ijk|V|lmn>"
           << endl;
      cout << i << "\t" << j << "\t" << k << "\t" << l << "\t" << m << "\t" << n
           << "\t" << J2ij << "\t" << J2lm << "\t" << J2 << "\t" << Tij << "\t"
           << Tlm << "\t" << T2 << "\t" << val << endl;

      cout << "a_p\t b_p\t c_p\t a\t b\t c\t J2ab_p\t J2ab\t T2ab_p\t T2ab\t "
              "J2\t T2\t "
           << endl;
      cout << a_p << "\t" << b_p << "\t" << c_p << "\t" << a << "\t" << b
           << "\t" << c << "\t" << J2ab_p << "\t" << J2ab << "\t" << T2ab_p
           << "\t" << T2ab << "\t" << J2 << "\t" << T2 << "\t " << val << endl;

      exit(0);
      continue;
    }

    double val_z = this->vec[ch_TP][tp_final][tp_init];

    out_file << i << "\t" << j << "\t" << k << "\t" << l << "\t" << m << "\t"
             << n << "\t" << J2ij << "\t" << J2lm << "\t" << J2 << "\t" << Tij
             << "\t" << Tlm << "\t" << T2 << "\t" << std::setw(8) << val
             << "\t";
    if (this->vec_flag[ch_TP][tp_final][tp_init] == 1) {
      out_file << std::setw(8) << val_z;
      double ratio = val / val_z;
      double diff = val - val_z;
      if (std::abs(ratio - 1) > 0.02 && std::abs(diff) > 0.00000000001)
        out_file << "\t * \t " << std::setw(3) << std::abs(ratio - 1) * 100
                 << "%";
    } else if (this->vec_flag[ch_TP][tp_init][tp_final] == 1) {
      val_z = this->vec[ch_TP][tp_init][tp_final];
      out_file << std::setw(8) << val_z;
      double ratio = val / val_z;
      double diff = val - val_z;
      if (std::abs(ratio - 1) > 0.02 && std::abs(diff) > 0.00000000001)
        out_file << "\t * \t " << std::setw(3) << std::abs(ratio - 1) * 100
                 << "%";
    } else {
      out_file << "null";
    }
    out_file << endl;
  }
  me_file.close();
  out_file.close();
}

int Package::FindTP_index(int a, int b, int c, int J2ab, int T2ab, int J2,
                          int T2) {
  int res = -1;
  int ch_TP = TP.FindC(J2, T2);
  int size = TP.state[ch_TP].size();
  for (int i = 0; i < size; i++) {
    if (a != TP.state[ch_TP][i].a || b != TP.state[ch_TP][i].b)
      continue;
    if (c != TP.state[ch_TP][i].c)
      continue;
    if (J2ab != TP.state[ch_TP][i].Jab_2 || T2ab != TP.state[ch_TP][i].Tab_2)
      continue;
    res = i;
    break;
  }
  return res;
}

void Package::print_TP_ALL(std::string filename, SingleP &SP, int h_core_Tz,
                           int p_beg_Tz, int p_end_Tz, double hw, int reg_pow,
                           int Lambda_Cut, vector<int> oc_o2z) {
  // y= oc_o2z[x] : y is state_Tz;  x is from others
  vector<int> oc_z2o;
  int size_oc = oc_o2z.size();
  oc_z2o.resize(size_oc);
  for (int x = 0; x < size_oc; x++) {
    int y = oc_o2z[x];
    oc_z2o[y] = x;
  }

  std::ofstream file;
  file.open(filename);
  LEC_CONST LECs;
  file << " ---- Parameters : ----" << endl;
  file << "\t Lambda_Cut = " << Lambda_Cut << "  (MeV)" << endl;
  file << "\t c1 = " << LECs.LEC_C1 << "  (GeV^-1)" << endl;
  file << "\t c3 = " << LECs.LEC_C3 << "  (GeV^-1)" << endl;
  file << "\t c4 = " << LECs.LEC_C4 << "  (GeV^-1)" << endl;
  file << "\t cD = " << LECs.LEC_CD << endl;
  file << "\t cE = " << LECs.LEC_CE << endl;
  file << "\t gA = " << LECs.LEC_gA << endl;
  file << "\t Hbar-Omega = " << hw << endl;
  file << "\t The Power of Regulator : " << reg_pow << endl;
  int h_core = std::floor(h_core_Tz / 2.0);
  int p_beg = std::floor(p_beg_Tz / 2.0);
  int p_end = std::floor(p_end_Tz / 2.0);
  /*
    file << "  === Single-Particle Orbits === " << endl;
    file << "\t index\t n\t l\t j2\t";
    for (int i = p_end; i >= p_beg; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Model Space" << endl;
    file << "\t --- --- --- --- \t";
    for (int i = h_core; i >= 0; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Core Space" << endl;
    file << endl; */

  file << "  === Single-Particle Orbits === " << endl;
  file << "\tindex\t n\t l\t 2*j\t";
  for (int i = p_end; i >= p_beg; i--) {
    int y = oc_o2z[i];
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state[y].n << "\t";
    file << SP.state[y].l << "\t";
    file << SP.state[y].j_2 << "\t";
    // file << SP.state[i].tz_2 << "\t";
  }
  file << "\t Model Space" << endl;
  file << "\t --- --- --- --- \t";
  for (int i = h_core; i >= 0; i--) {
    int y = oc_o2z[i];
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state[y].n << "\t";
    file << SP.state[y].l << "\t";
    file << SP.state[y].j_2 << "\t";
    // file << SP.state_Tz[i].tz_2 << "\t";
  }
  file << "\t Core Space" << endl;
  file << endl;
  file << "a\tb\tc\td\te\t f\t J2ab\t J2cd\t J2\t T2ab\t T2cd\t T2\t\t "
          "<abc|V|def>"
       << endl;
  int ch_size = TP.channel.size();
  for (int TP_ch = 0; TP_ch < ch_size; TP_ch++) {
    int size = TP.state[TP_ch].size();
    for (int i = 0; i < size; i++) {
      int a = TP.state[TP_ch][i].a;
      int b = TP.state[TP_ch][i].b;
      int c = TP.state[TP_ch][i].c;
      int J2ab = TP.state[TP_ch][i].Jab_2;
      int T2ab = TP.state[TP_ch][i].Tab_2;
      int J2 = TP.state[TP_ch][i].J_2;
      int T2 = TP.state[TP_ch][i].T_2;
      for (int j = 0; j < size; j++) {
        int d = TP.state[TP_ch][j].a;
        int e = TP.state[TP_ch][j].b;
        int f = TP.state[TP_ch][j].c;
        int J2de = TP.state[TP_ch][j].Jab_2;
        int T2de = TP.state[TP_ch][j].Tab_2;
        double val = vec[TP_ch][i][j];
        if (std::abs(val) > 0.0000000000001) {
          int a_x = oc_o2z[a];
          int b_x = oc_o2z[b];
          int c_x = oc_o2z[c];
          int d_x = oc_o2z[d];
          int e_x = oc_o2z[e];
          int f_x = oc_o2z[f];
          file << a << "\t " << b << "\t " << c << "\t ";
          file << d << "\t " << e << "\t " << f << "\t ";
          file << J2ab << "\t " << J2de << "\t " << J2 << "\t ";
          file << T2ab << "\t " << T2de << "\t " << T2 << "\t\t ";
          file << val << endl;
          file << endl;
        }
      }
    }
  }
  file.close();
}

void Package::print_TP_ALL_f(std::ofstream &file, SingleP &SP) {
  // std::ofstream file;
  // file.open(filename, ios::std::app);
  // LEC_CONST LECs;
  // file << " ---- Parameters : ----" << endl;
  // file << "\t Lambda_Cut = " << Lambda_Cut << "  (MeV)" << endl;
  // file << "\t c1 = " << LECs.LEC_C1 << "  (GeV^-1)" << endl;
  // file << "\t c3 = " << LECs.LEC_C3 << "  (GeV^-1)" << endl;
  // file << "\t c4 = " << LECs.LEC_C4 << "  (GeV^-1)" << endl;
  // file << "\t cD = " << LECs.LEC_CD << endl;
  // file << "\t cE = " << LECs.LEC_CE << endl;
  // file << "\t gA = " << LECs.LEC_gA << endl;
  // file << "\t Hbar-Omega = " << hw << endl;
  // file << "\t The Power of Regulator : " << reg_pow << endl;
  //
  // file << " ---- 3BME range ----  " << endl;
  // int num_t = Same_particle_Flag_g.size();
  // for (int i = 0; i < num_t; i++) {
  //   int Same_particle_Flag = Same_particle_Flag_g[i];
  //   if (Same_particle_Flag > 0) {
  //     file << "!! The index of <a'b'c'|V|abc> has the index constrain !!"
  //          << endl;
  //     file << "    ";
  //     if (Same_particle_Flag == 1)
  //       file << "\t a' = a " << endl;
  //     if (Same_particle_Flag == 2)
  //       file << "\t b' = b " << endl;
  //     if (Same_particle_Flag == 3)
  //       file << "\t c' = c " << endl;
  //     if (Same_particle_Flag == 12)
  //       file << "\t a' = a & b' = b" << endl;
  //     if (Same_particle_Flag == 13)
  //       file << "\t a' = a & c' = c" << endl;
  //     if (Same_particle_Flag == 23)
  //       file << "\t b' = b & c' = c" << endl;
  //     if (Same_particle_Flag == 123)
  //       file << "\t a' = a & b' = b & c' = c" << endl;
  //   }
  //   SP.print(TP_allow_vec_g[i], file);
  // }

  double EPS = 1E-8;
  file << " Only output abs(ME) > " << EPS << endl;
  file << "a'  b'  c'  a  b  c   J2ab  J2cd  J2  T2ab  T2cd  T2 \t "
          "<abc|V|def>"
       << endl;
  int ch_size = TP.channel.size();
  for (int TP_ch = 0; TP_ch < ch_size; TP_ch++) {
    int size = TP.state[TP_ch].size();
    for (int i = 0; i < size; i++) {
      int a = TP.state[TP_ch][i].a;
      int b = TP.state[TP_ch][i].b;
      int c = TP.state[TP_ch][i].c;
      int J2ab = TP.state[TP_ch][i].Jab_2;
      int T2ab = TP.state[TP_ch][i].Tab_2;
      int J2 = TP.state[TP_ch][i].J_2;
      int T2 = TP.state[TP_ch][i].T_2;
      for (int j = 0; j < size; j++) {
        int d = TP.state[TP_ch][j].a;
        int e = TP.state[TP_ch][j].b;
        int f = TP.state[TP_ch][j].c;
        int J2de = TP.state[TP_ch][j].Jab_2;
        int T2de = TP.state[TP_ch][j].Tab_2;
        // double val = vec[TP_ch][i][j];
        double val = vec_S.val(TP_ch, i, j);
        if (std::abs(val) > EPS) {
          file << a << "  " << b << "  " << c << "  ";
          file << d << "  " << e << "  " << f << "  ";
          file << J2ab << "  " << J2de << "  " << J2 << "  ";
          file << T2ab << "  " << T2de << "  " << T2 << "\t ";
          file << val << endl;
          // file << endl;
        }
      }
    }
  }
}
